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Abstract. By numerical modeling we investigate fluid transport in low-Reynolds-number flow achieved 
with a special elastic filament or artifical cilium attached to a planar surface. The filament is made of 
superparamagnetic particles linked together by DNA double strands. An external magnetic field induces 
dipolar interactions between the beads of the filament which provides a convenient way of actuating the 
cilium in a well-controlled manner. The filament has recently been used to successfully construct the first 
artificial micro-swimmer [R. Dreyfus at ai, Nature 437, 862 (2005)]. In our numerical study we introduce a 
measure, which we call pumping performance, to quantify the fiuid transport induced by the magnetically 
actuated cilium and identify an optimum stroke pattern of the filament. It consists of a slow transport 
stroke and a fast recovery stroke. Our detailed parameter study also reveals that for sufficiently large 
magnetic fields the artificial cilium is mainly governed by the Mason number that compares frictional to 
magnetic forces. Initial studies on multi-cilia systems show that the pumping performance is very sensitive 
to the imposed phase lag between neighboring cilia, i.e., to the details of the initiated metachronal wave. 

PACS. 87.19.St - 87.16.Ac - 87.16.Qp 



1 Introduction 

Fluid transport and mixing on the microscopic level at 
J low Reynolds numbers is a fascinating problem that is at 

■ the center of a successful lab-on-chip technology [I]. Na- 
' ture has provided an ingenious solution to this challenge 
. by using long elastic filaments, called flagella or cilia, that 

■ are actuated internally by molecular motors '2,3,4]. The 
' resulting beating patterns, in general three-dimensional, 
' have to be non-reciprocal as Purcell has taught us i5j . Na- 

■ ture uses arrays of collectively beating cilia to transport 
mucus in the respiratory tract, fluid in the brain 0, or to 
propel microorganisms such as the Paramecium. During 
an early stage of a developing embryo, arrays of rotating 
cilia are responsible for establishing the left-right asymme- 
try in the placement of organs [7]. Recently, experimen- 
tal efforts have been initiated to copy nature's successful 
concept by developing biomimetic or artificial cilia that 
are actuated by external fields [DIllll] or to move fluid 
with the help of bacterial carpets [10] . On the other hand, 
there has been an increasing interest in recent times in 
contributing to the theoretical understanding of how sin- 
gle cilia or flagella function (e.g. Refs. [TT1IT2l[T3l[T4l[T5l[T6l 
[TTlfTMro] ) and of how their collective beating patterns, 
known as metachronal waves, occur (e.g. Refs. p!2l[T7ll20l 
[2T1[22l[28l[24li25.2fiJ). 

Some years ago, Dreyfus et al. introduced the first 
artificial micro-swimmer [8] based on an artificial flagel- 



lum. This elastic filament consists of superparamagnetic 
micron-sized beads that are linked together by pieces of 
double-stranded DNA p7l[28]. The flagellum is actuated 
by an oscillating magnetic field so that it drags an at- 
tached red-blood cell forward. Modeling using a contin- 
uum [T5] or a discretized [18j approach described the be- 
havior of the swimmer very well. 

Based on our theoretical work in Ref . [18] , we explore 
in this article by numerical modeling how the artificial 
fiagellum or cilium can be employed for fluid transport. We 
attach the cilium to a surface and actuate it by an external 
magnetic field, the direction of which oscillates about the 
surface normal in an asymmetric fashion. To achieve fluid 
transport, the stroke pattern has to be asymmetric |17j . 
We therefore introduce a slow transport stroke, where the 
cilium remains nearly straight, and a fast recovery stroke, 
where the cilium bends due to increased hydrodynamic 
friction. A measure for the amount of transported fluid 
during one beating cycle, called pumping performance, 
helps us to discuss the system in detail and to identify 
the optimum conditions under which the cilium should be 
operated. We also look at multi-cilia systems with a de- 
fined phase lag between neighboring cilia and illustrate 
that such a phase lag is advantageous for fluid trans- 
port. Of course, our metachronal waves are controlled by 
the external field and do not occur through self-organized 
synchronization of the beating cilia. There is now strong 
evidence that this intriguing feature of biological cilia ar- 
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Fig. 1. Filament attached perpendicular to a bounding wall. 
The first three bonds to and ti are shown. The actual 
filament starts with the bead at ro, whereas the virtual beads 
drawn as dotted lines are kept at fixed positions to anchor the 
filament perpendicular to the xy plane through their elastic 
force contributions. The actuating magnetic field oscillating 
around the z axis is also sketched. 

rays is caused by hydrodynamic interactions [T2l[T7l[20ll2Tl 
[^2112511^ Il251l26j ) . Similar synchronization phenomena me- 
diated by hydrodynamic interactions were discussed for 
sperm cells by Taylor [3D] and for helical flagella in Ref. 

The article is organized as follows. Section [2] explains 
our modeling of the artificial cilium including the treat- 
ment of hydrodynamic interactions close to surfaces and 
the magnetic actuation cycle of the filament. Section [3] in- 
troduces the pumping performance as a measure for the 
transported fluid. In Sections 2] and [5] we then discuss the 
pumping performance for a single cilium and for multi- 
cilia systems. The article ends with a conclusion. 



motion. Details of the derivations are given in Ref. |18j . 
The equations of motion contain hydrodynamic mobilities 
which we construct up to Rotne-Prager level based on the 
appropriate Green function commonly called Blake's ten- 
sor [551fM] . This guarantees a vanishing fluid velocity field 
at the bounding surface. The description of the time pro- 
tocol of the actuating magnetic field concludes the section. 



2.1 Energies and forces of the superparamagnetic 
filament 

The total free energy for our bead-spring model is given by 
a sum over dipolar and elastic-energy contributions [18] . 

H = + + , (1) 

where is the dipole-dipole interaction energy, is 
the stretching energy obeying Hooke's law and is the 
discretized free bending energy of an elastic rod. As in our 
previous work, the actuation and hence also the motion of 
the filament takes place in the yz plane and no twisting of 
the filament occurs, allowing us to ignore a twisting energy 
in our model. We shall briefly discuss the individual energy 
contributions of Eq. ([T]) in the following. 

Let ti be the vector connecting the centers of the two 
adjacent beads labeled i and i — A deviation k = \ti \ of 
the beads from their equilibrium spacing Iq then gives a 
total stretching free energy of 

1 

2 = 

where k is the stretching constant and N the total number 
of beads. 

Discretizing the continuum bending energy of the worm- 
like chain model [5^1155] . we obtain by a straightforward 
calculation [TB] 



2 Model 

The superparamagnetic fllament is modeled by a bead- 
spring conflguration, which additionally resists bending 
like a worm-like chain [32]. Consequently, each bead of 
the filament is subject to stretching and bending forces, 
for which the chemical linkers are responsible, and to dipo- 
lar interaction forces due to the induced magnetic dipoles 
of the beads. We completely ignore the contributions of 
the chemical linkers to hydrodynamic friction, so the fil- 
ament interacts with the fiuid surrounding solely via the 
hydrodynamic friction of the beads. As illustrated in Fig. 
[TJ the filament is attached to a planar surface with the 
help of two virtual beads that fix its position and give 
it an orientation orthogonal to the surface. These virtual 
beads contribute to elastic forces but do not participate 
in hydrodynamic or dipolar interactions. 

In the following, we will summarize the forces, acting 
on each bead within the filament, and the equations of 



A 

'0 , 



H+l 



(3) 



where ii ~ ti/li and A is the bending stiffness. 

Finally, consider two identical paramagnetic beads with 
radius a and magnetic susceptibility x subject to a homo- 
geneous external magnetic field B. Both beads develop a 
dipole moment with identical orientation and strength. 



P 



3^0 



XB. 



(4) 



where fiQ = 4tt x lO^^N/A^ is the permeability of free 
space. The resulting dipoles of the beads of the filament 
give rise to the total dipole-dipole interaction energy 



3(p-ry) 



9^0 



(5) 
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where = \rj — nj, = {rj — ri)/rij and the double 
sum runs over all terms with i < j, where the minimal i 
is and the maximal j is iV — 1. 

Now that all relevant interaction energies have been 
characterized, the non-hydrodynamic force acting on each 
bead is obtained as usual, 



(6) 



where Vr-i is the gradient operator with respect to r^. The 
total non-hydrodynamic force Fi = + F^ + F-^ acting 
on bead i is readily calculated from Eqs. (I2l[3l[5]). Explicit 
expressions of Fi can also be found in Ref. ^ISj . 



2.2 Equations of Motion 

On the micron length scale, the motion of a particle im- 
mersed in a viscous fluid, such as water, is entirely domi- 
nated by friction and the particle's inertia can be neglected 
provided the time scale of interest exceeds the momentum 
relaxation time Hence, the velocities Vi of the beads 
are proportional to the forces Fj acting on them and the 
beads obey the following equations of motion [36] 



■Wi = ^ /Xy with Fj = 



(7) 



All the forces depend on the spatial configuration of the 
filament, i. e., the beads' locations r; . Furthermore, the 
dipolar forces also possess an explicit time dependence 
through the external magnetic field, which we use to ac- 
tuate the filament. 

Hydrodynamic friction enters the equations of motion 
via the mobilities /Xy, which depend on the geometrical 
configuration of the beads. The flow field induced by one 
moving bead creates a drag on all other beads in the 
vicinity and thus indirectly displaces them. Since induced 
flow fields are long ranged (they decay as l/r, where r 
is the distance from a moving bead), hydrodynamic in- 
teractions play an important role in viscous systems with 
low Reynolds number. While hydrodynamic interactions 
constitute a highly complicated many-body problem [36] , 
their leading order is given by two-particle interactions. 
The relatively simple Rotne-Prager approximation can be 
employed whenever the beads are not so close together 
that lubrication effects need to be considered [36. 37^ . How- 
ever, the surface, to which our filament is attached, in- 
troduces significantly more complexity due to the no-slip 
boundary condition even when we restrict ourselves to 
pairwise interactions at the Rotne-Prager level. 



2.3 Hydrodynamic interactions 

In the low Reynolds number regime, the fluid flow velocity 
at an arbitrary point r is linearly related to a point force 
Fq at r' by |3S] 



, , Gir- r') ^ 



(8) 



where u{r) denotes the flow field at r and r] is the fluid's 
viscosity. In an unbounded fluid, the Green function G{r~ 
r'), commonly referred to as Oseen tensor in literature, is 
given by 



G(r) 



1 



(9) 



Here, / is the 3x3 identity matrix and (gD denotes the 
dyadic product. 

The aforementioned Rotne-Prager mobilities of spheri- 
cal particles are commonly used quantities in an unbounded 
fluid, which can be derived from the Oseen tensor [36]. 
The respective self and cross mobilities are given by the 
following expressions for two spheres of radius a, 



(10) 



, ^^J , (11) 



where 



rj and rij 



j/vij and /io = (671770) ~ 



However, in close proximity to a planar surface with 
no-slip boundary condition, the traditional Rotne-Prager 
mobilities can no longer be employed. The pressure and 
velocity fields of a point force for this boundary condition 
have been known for a long time and were first derived by 
Lorentz [40]. Blake put these results into a modern form 
replacing the Oseen tensor by the appropriate Green func- 
tion |33J. The condition of a vanishing fluid velocity field 
on an infinitely extended plane is satisfied with the help of 
appropriate mirror images, similar to the image charge ap- 
proach often used in electrostatics. In contrast to electro- 
statics, where it suffices to simply mirror the charge distri- 
bution, the hydrodynamic image system is more involved 
due to the more complicated structure of the Stokes equa- 
tion when compared to the Poisson equation. Hence, so- 
called stresslet and source-dipole contributions are needed 
in addition to the stokeslet of the point force and its mir- 
rored point disturbance (also called anti-stokeslet). This 
yields Blake's tensor, 

= G{r^r') 

-G{r-r') + 5G'"'{r,r') , (12) 

where G{r — r') and G(r — r') are Oseen tensors, r' is the 
coordinate vector of the stokeslet source, and r' is the posi- 
tion of the anti-stokeslet source, i. e., the stokeslet source 
mirrored at the bounding xy plane [r — {x,y,z),r' = 
(cc', y', —z')]. Finally, 5G™(r, r') denotes the source-dipole 
and stresslet contributions. 

For particles that are far apart from each other, it suf- 
fices to use the Blake tensor for the mobility functions, 
i. e., to treat the particles as point-like objects, as is fre- 
quently done [23,41J. However, if the particles approach 
each other, their finite sizes become relevant as is the case 
in our filament. To take this effect into account, we derive 
the mobility functions in Eq. ([7|) up to the Rotne-Prager 
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level following Refs. [42jl43j. Consider f{r') to be the force 
density on the surface dVj of a sphere j with its centre at 
rj . The force density acts on the fluid and induces the flow 
fleld 

uir) = ^[ G^'^^^(r,r')/(r')d5' (13) 

The leading order of f{r') is Fj/ (47ra^), where Fj denotes 
the total force on particle j. Taking it into account and 
expanding QBiake^j,^^,/-^ about r' = rj yields for the flow 
field 



u{r) « 



6 



(14) 

We proceed by calculating the effect of u{r) on the motion 
of another sphere i with the help of Faxen's theorem |,36 , 



F, 



1 



u{ri) 



(15) 



Note that u(ri) also includes the flow fields created by 
the mirror image of sphere i, i.e., by the term G'^™'{r,r') 
in the Blake tensor of Eq. (fT2|) . Combining Eqs. (fT4|) and 
(fT5|) . we immediately obtain the cross mobilities ^44] 



1 



1 



1 



and the self mobility of particle i, 



G^''^'=^(n,r,) 
(16) 



Mr 



1 

H 



±1 

3a 



6 



1 H V 

6 ^ 



G™(r„r,) 



(17) 



Application of the differential operators in Eqs. and 
(fT7|) to the Oseen tensors G{r — r') and G{r — r') in the 
Blake tensor simply gives Rotne-Prager matrices. Hence, 
we write the mobilities as 



Mm: 



- M'^^ln - Ti) + SUself 



(18) 
(19) 



where /x''P(r) are the Rotne-Prager tensors defined by Eq. 
(fTT|) and Sfiseif and J/Xy are the sourcelet and stresslet 
contributions originating from SG'^"^. These evaluate to 
more complicated expressions which we give in Appendix 



2.4 Reduced Equations of Motion 

We have shown in previous work that the number of sys- 
tem parameters for the superparamagnetic filament can 
be significantly reduced by introducing just three dimen- 
sionless variables . The essential parameters governing 
the dynamics of the filament were identified by rescaling 
the dynamic equations appropriately so that only reduced 
variables appeared. Instead of repeating the derivation of 



these reduced variables given in Ref. [18j . we limit our- 
selves to discussing their physical significance. 
One of the emerging quantities is 



A 



1/4 



L 



(20) 



which we call 'sperm number' in the following [14j . The 
physical interpretation of Sp is that it compares bend- 
ing to frictional forces in our bead-spring chain. It bears 
a resemblance to a dimensionless quantity well-known 



from analytical continuum models [TT|fT4j when substitut- 
ing Gnrja/lo by 7x which is the perpendicular friction con- 
stant per unit length of a slender body. Note that dirria/lo 
approximates 7^ for our bead-spring chain when hydrody- 
namic interactions are neglected. Via Sp = L/lh, one can 
assign a characteristic length to the system. In the con- 
tinuum model of a slender body, the analogous quantity 
is the elastohydrodynamic penetration length [TT] . An 
oscillation with frequency lo started at one of the ends of 
a sufficiently long filament penetrates the filament up to 
the length l^. Conversely, if the filament is much shorter 
than Z^, i. e., if L ^ l^^ the filament oscillates like a rigid 
rod over its whole length. 

A second important quantity is the reduced magnetic 
field strength. 



= 779 ~^ ■ 



(21) 



which more suitably describes the effect of the external 
field on the superparamagnetic filament. The number 
compares dipolar to bending forces and it is proportional 
to the magnetoelastic number introduced in Ref. [51IT5]. 
Finally, we obtain a reduced spring constant 



(22) 



While Sp and Bs are experimentally accessible via the 
frequency uj and the magnetic field strength _B, fcg is a 
fixed quantity of the superparamagnetic filament itself. 
Therefore, we do not make explicit use of the stretching 
mode of the filament and we always choose a sufficiently 
large kg to keep overall length fiuctuations well below 10%. 

In what follows, we find that in the interesting regime 
for fluid transport the magnetic forces of the induced dipoles 
dominate over the bending forces of the superparamag- 
netic filament. It therefore makes sense to introduce a 
fourth dimensionless number as the ratio of frictional to 
magnetic forces. In literature on magnetorheological sus- 
pensions, it is also called the Mason number [45.46J: 



Sp/Bl 



(23) 



2.5 Actuation of the filament 



The direction of the actuating magnetic field of strength 
B oscillates in the yz plane (see Fig. [21): 



B{t) = {Q,B sin ip{t),B cos (p{t)) , 



(24) 
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Fig. 2. The angle y enclosed by the magnetic field B{t) and 
the z axis is shown as a function of time, ip has different 
velocities when decreasing and increasing. The relative dif- 
ference between and r; defines the asymmetry parameter 
e = {ti - Ts)/{ts + n). 



where if{t) is the angle the field encloses with the z axis. 
In our modeling, the angular amplitude ipmax was always 
60°. The magnetic field induces a motion of the filament 
caused by the time-dependent dipole-dipole interaction of 
the filament's constituent beads. 

In order to accomplish net fluid transport, a non-rc- 
versible motion of the filament is required [5] . In addition, 
the motion of the filament in positive and negative y direc- 
tion has to be asymmetric to achieve fluid transport along 
the y axis. In nature, the beating cycle of biological cilia 
can typically be broken down into a power or transport 
stroke, during which the cilium moves in a fashion similar 
to a rigid rod with high hydrodynamic resistance, and a 
recovery stroke, where it bends strongly and moves back 
with reduced resistance close to the surface. We mimic 
this behavior by introducing a faster stroke in one direc- 
tion followed by a slower stroke back to the original po- 
sition to complement the beating cycle [47j. However, in 
contrast to real cilia the slower stroke of our fllament is 
the transport ('power') stroke and the faster stroke serves 
as the 'recovery' with a bent fllament. We will illustrate 
this below. 

To accomplish the desired asymmetry, we simply ro- 
tate the B vector fast from one side to the other and move 
it back more slowly. Let T = 2Tr/u! be the duration of a 
full actuation cycle. As illustrated in Fig. [21 we split T into 
two parts: T = Ts + Ti. During the fast stroke with dura- 
tion Ts, If is decreasing with velocity —2ipmax/Ts and ip 
is increasing during the slow stroke with duration r/ with 
velocity 2ipmax/Ti- To quantify the asymmetry in the ac- 
tuation cycle and therefore in the beating pattern of the 
fllament, we define the asymmetry parameter 



Tl 



(25) 



which is zero for Ts = ti and tends to one in the limit of 

Tl > Ts- 

In comparison with biological cilia or the work of Kim 
et al. [T7], the magnetic actuation technique has the ad- 
vantage of not requiring any active elements within the 



filament or its base, such as a driving motor with a geo- 
metrical switch for force reversal. Indeed, it entirely suf- 
fices to manipulate the more easily accessible macroscopic 
control field, whose time protocol allows us to achieve an 
asymmetric and non-reciprocal beating cycle. 



3 Measuring the pumping performance 

While a non-reciprocal and asymmetric motion of the flla- 
ment is required to generate a net fluid transport, a quan- 
tity to measure the pumping performance of the fllament 
is far from obvious. In the following we will introduce such 
a measure based on the works in Refs. PTIHS] . 

The main idea for the pumping performance is to inte- 
grate the fluid flow initiated by the beating fllament over 
a whole plane parallel to the bounding surface. The inte- 
grated fluid flow is determined by the laterally averaged 
Blake tensor which assumes a particularly simple form 
HHl: 



dxdyGfy''''^{x,y,z,z') 



z -\- z' — \z — z'\ min(2;, z') 



21] 



V 



(26) 



where we restricted ourselves to the yy component since 
fluid is transported along the y direction. Then, the inte- 
grated flow generated by all the beads of the fllament 
is approximated by }49| 



N 



' 1— n 



(27) 



where Zi is the distance between the wall and bead i and 
Fyi is the force acting on it in y-direction. This expression 
is a function of the distance z from the wall and it is con- 
stant for z > Zi{i = 1...N). For our beating filament, J-'{z) 
changes with time but is cyclic over a full beating cycle 
of the filament as is illustrated in Fig. [3K) for z > L. Fig- 
ure [31d) shows !F{z) as a function of z at different points 
of time during a typical beating cycle. As expected, J-{z) 
only depends on z roughly over the length L of the flla- 
ment and then becomes constant. Now, we consider the 
time average of T{z > L) over one beating cycle as a suit- 
able measure for fluid transport, which we call pumping 
performance: 



^ -1 



t+T 



dt' T{z > L,t') , 



(28) 



where T is the period of one actuation cycle. 

To understand how well the magnetically actuated flla- 
ment transports fluid, we compare it to an idealized stroke 
of a rigid rod with the same length parameter L = {N — 
l)lo and the same thickness as the fllament (see Fig.|4]): 

1. In the transport stroke the rod is oriented perpendic- 
ular to the bounding surface and it is dragged parallel 
to the surface along a distance L keeping its center a 
distance L/2 + 3a above the surface. 
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1.5 ^ 



2.5 3 3.5 
time (in units ofr ) 



-6 -4-2 2 
integ. {lowJ-(z) (a.u.) 




0.4 0.6 
d (in units of L ) 

Fig. 3. a) Integrated flow !F{z > L) in arbitrary units as a 
function of time in units of T for the parameters Sp = 3, Bs = 3 
and e — 0.9. b) Integrated flow T{z — 1.5L) as a function of 
height z at different points of time during one beating cycle. 
The stroboscopic snapshots to the right correspond to the slow 
part of the stroke and the ones to the left to the faster stroke, c) 
Pumping performance (0) times period T of the idealised 
reference cycle as a function of the parameter d as deflned by 
Eq. (|29p . The y axis is scaled to units of the maximal value 
(0) = 3.895 X 10"" m^/s at d = 0. 



1) transport stroke 



2) recovery stroke 



CXXKX3 



, lo = 3a 



Fig. 4. Geometry of the idealized stroke of a rigid rod con- 
sisting of A'' particles of radius a with distance lo between the 
centers of adjacent beads. The length parameter L is deflned 
as L = (iV - l)lo 



2. The rod is then rotated by 90° to be parallel to the sur- 
face and then in the recovery stroke it is dragged along 
its long axis to its original position with a separation 
of d + 3a above the surface. 

In unbounded Stokes flow, the first and second part of the 
stroke use the respective maximal and minimal hydrody- 
namic resistance to produce a net fluid transport over one 
cycle. In the presence of a no-slip wall, the hydrodynamic 
resistance increases closer to the wall. Therefore, the dis- 
tance d in the second part turns out to be an important 
parameter. 



For the idealized stroke, we define a reference pumping 
performance 



t+T 



dt' T{z > L,t') , 



(29) 



where T — T± + Tn is the time required for both parts 
of the cycle and the variable d emphasizes the separation 
between the rigid rod and the surface for the second part 
of the reference cycle. We determined the reference pump- 
ing performance with the help of our filament as follows. 
We let a constant force act on each bead and numerically 
evaluated the average velocity of the beads in both parts 
of the strokes from which we then calculated the times T± 
and Til the filament takes to move the length L. Together 
with Eq. ([27| T^^{d) is then determined. In particular, 
it becomes clear that !FJ^^{d)T only depends on effec- 
tive friction coeffcients, L, and d. Note that in both parts 
of the reference stroke, the filament experiences an addi- 
tional torque which alters its orientation. To really realize 
the idealized stroke, one would have to counterbalance the 
frictional torques appropriately. Fig[3t) shows TJ^f{d)T 
as a function of d. A linear fit nicely connects the data 
points acquired by simulations and represents the linear 
dependence of !F{z > L,t') on d in the recovery stroke. 
It makes sense that maximal pumping occurs when the 
recovery stroke of the reference cycle is at a minimum dis- 
tance 3a from the wall [d = 0) because the no-slip bound- 
ary condition reduces the fluid volume the rod can drag 
along. Interestingly, a reversal of the pumping direction 
sets in at about d — 0.9. Here, the perpendicular move- 
ment of the rod gives rise to less fluid transport compared 
to the parallel one simply because it is on average much 
closer to the boundary. This highlights the importance of 
hydrodynamic interactions not only amongst the beads 
but also with the bounding wall. 

The reference stroke describes an optimum fluid trans- 
port one can achieve with a rigid rod of length L close to a 
surface. We therefore compare the pumping performance 
J-oo of the magnetically actuated filament to the reference 
value T^^^ [d = 0) at d = i.e., when the ideahzed refer- 
ence cycle assumes ist maximum value. This results in the 
reduced pumping performance 



(30) 



^-/(d = 0) 

which we use in the following. Expression Too in Eq- (|30p 
is evaluated after the filament has assumed its limiting cy- 
cle, which typically happens after less than five actuation 
cycles after starting the simulation. We have checked that 
^ evaluates to zero either if the motion of the filament is 
symmetric (r^ — ti) or if it is reciprocal, i.e., if the fila- 
ment does not bend, e.g., in the case Sp — > 0. Hence, Eq. 



([SO]) constitutes a suitable measure for characterizing the 
pumping performance of the magnetically actuated artifi- 
cial cilium. 
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Fig. 5. Reduced pumping performance ^ (a) and the time- 
averaged dissipated power v in arbitrary units (b) for a single 
cilium as a function of sperm number and asymmetry pa- 
rameter e. The reduced magnetic field strength is Bs ~ 2.5. 
The white dots mark parameters for which stroboscopic snap- 
shots of the filament are shown in Fig. 

4 Pumping performance of a single artificial 
cilium 

Figure [5^) shows the pumping performance ^ of a single 
filament as a function of the Sperm number Sp and the 
asymmetry parameter e at a fixed strength of the magnetic 
field Bs = 2.5. The most striking feature is that the perfor- 
mance is strongly peaked for e close to one and ai Sp ~ 3. 
Such a peak is also observed in the swimming velocity of 
the artificial micro-swimmer [51 ll51fT5] . The corresponding 
stroke pattern for Sp ~ S is illustrated in the middle pic- 
ture of Fig. IHK) . In the slow transport stroke the filament 
rotates clockwise being nearly straight. It uses the high 
friction coefhcient of a rigid rod dragged perpendicular to 
its axis to pump fluid. In the fast recovery stroke, the fila- 
ment bends due to the large hydrodynamic friction forces 
that scale with velocity and then relaxes back to the ini- 
tial configuration. As Fig. [5^) illustrates, fluid transport 
is also noticeable in the recovery stroke. So, the pumping 
performance, even for the most efficient stroke pattern, is 
the result of a small asymmetry in the amount of fluid 
transported to the right or left. In the example of Fig. 



[3K), which is close to the optimum stroke pattern, only 
4.3 % of the total amount of moved fluid are effectively 
transported in positive y direction. As a result, the maxi- 
mum pumping performance in Fig. [5^) is only 6% of the 
reference stroke. As expected, the pumping performance 
vanishes for symmetric beating about the z axis, i.e., when 
e = 0. The same is true for Sp Q: The filament follows 
the actuating magnetic field instantaneously. Aside from 
the base, the filament therefore remains straight and the 
stroke is reciprocal. The left picture in Fig. [6^) gives an 
example of such a stroke. A reversal of the pumping direc- 
tion < 0) occurs at Sp ~ 5.5, albeit only with a rather 
weak performance. Finally, the pumping performance goes 
to zero for increasing Sp or frequency since the filament 
can no longer follow the actuating field as illustrated in the 
right picture of Fig. [5^). In summary, an optimal pumping 
performances is only achieved for intermediate values of 
Sp. 

Let us add a further remark. Naively, one might antic- 
ipate that the closer e is to one, the more pronounced the 
incurred bending during the fast part of the actuation cy- 
cle, giving rise to better pumping performance. However, 
a limit on the maximal speed with which the filament 
manages to follow the field exists. The extreme case of 
instantaneous switching of the direction of the magnetic 
field in the recovery stroke (r^ = 0) is not conducive to 
generating fiuid transport. On the contrary, it seems im- 
portant that the free end of the filament follows the mag- 
netic field's orientation during the fast stroke. Relaxation 
from the bent to a straight configuration then happens 
during the time interval t; allotted to the slow part of the 
cycle. Consequently, the filament takes significantly longer 
to perform the actual recovery stroke of the cycle than the 
magnetic field. Naturally, this observation translates into 
the fact that the e value at which maximum pumping per- 
formance occurs decreases with increasing Sp. In Fig. [5] a), 
where e ranges up to 0.94 this is visible to the right of the 
optimal pumping performance. 

Part b) of Fig. O shows the time average of the dissi- 
pated power, 

1 rt+T N 

^ = ^/ E^'--^^^*'' (31) 

as a function of Sp and e. Surprisingly, it does not dis- 
play such a pronounced behavior as the pumping perfor- 
mance. Particularly, there is hardly any dependence on the 
asymmetry parameter e visible. Being proportional to the 
square of the beads' velocities, one could assume 
and therefore v ^ Sp (recall that Sp ^ uj^^^). At small 
Sp, the data shown in Fig. \Ejp) do show a steep incline. 
However, they do not increase as S^, because already at 
Sp around 2 the filament does not completely follow the 
actuating field as the snaphots in Fig. [6^) illustrate. This 
effect becomes more pronounced for increasing Sp, and the 
data clearly deviate from the naive assumption of u Sp. 
The third picture in Fig. [S^) demonstrates that already at 
Sp — 5 the filament lags significantly behind the actuating 
magnetic field. 
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a) fixed 8^=2.5 




b) fixed Sp =3 




Bs=l Bs^3 Bs=6 



Fig. 6. Stroboscopic snapshots of the filament at diff'erent 
times during the beating cycle for e = 0.9. The trajectory of 
the top bead during one beating cycle is also indicated. In the 
slow transport stroke the filament rotates clockwise, the fast 
recovery stroke occurs to the left, as indicated by the arrows. 
A pronounced bending of the filament occurs only at interme- 
diate sperm number Sp and magnetic field strength Bg. 



As illustrated by the three-dimensional plot of Fig. 
[TK) , there exists a pronounced dependence of the pumping 
performance ^ on the strength Bg of the actuating mag- 
netic field which has also been observed in the swimming 
velocity of the artificial micro-swimmer [8|I151IT8]. Figure 
[T}^) shows this behavior for different constant values of 
Sp. Increasing the magnetic field from zero, the pump- 
ing perfomance stays close to zero and then, beyond a 
certain threshold value, it grows until it reaches a maxi- 
mum. It finally decreases and even becomes negative. The 
snapshots in Fig. [6h) again explain this behavior. Small 
field strengths Bg are not high enough to overcome the 
hydrodynamic friction forces and therefore the motion of 
the filament is very limited. On the other hand, at large 
strengths B^ the filament is always straight and therefore 
performs a reciprocal motion. An optimal stroke only ex- 
ists in an intermediate regime for the strength Bs ■ Clearly, 
the optimal performance shifts with increasing Sp to larger 
values of Bg since a larger field is needed to move the fil- 
ament through the fluid. In other words, for larger Bs 
the filament is stiffer and a higher frequency to cx[ Sp is 
needed to achieve the optimum stroke. Higher frequency 
means larger frictional forces and therefore a larger pump- 
ing performance, as Figs. [TK) and b) demonstrate. When 
the magnetic forces on the filament exceed the bending 
forces, one expects the dynamics of the filament to be 
determined by the ratio of the hydrodynamic friction to 
magnetic forces, which we introduced in Eq. as Ma- 
son number Ma- From Fig. [7^) we extract curves £,[Sp) 
for different constant and rescale each of these curves 
with the respective maximum values In Fig. [Tt), all 
the data points for Bg > 2 are plotted as a function of the 
Mason number Ma- Indeed, most of the data points fall 
on a master curve, as predicted. Deviations occur for data 
points with Bg close to 2. 
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Fig. 7. a) Pumping performance ^ as a function of Sp and Bs 
for £ = 0.9. b) Plot of ^ versus Bs for different 5*^. The black 
dots mark parameters for which stroboscopic snapshots of the 
filament are shown in Fig.lBJs). c) Pumping performance ^ as a 
function of the Mason number Ma — Sp/Bs, where ^ is given 
in units of the maximum values when Bs is kept constant. 

5 Pumping performance of hydrodynamically 
interacting cilia 

Nature often uses arrays of beating cilia rather then a 
single isolated cilium for generating fluid transport or to 
propel microorganisms such as a Paramecium [2]. These 
cilia beat in a synchronized fashion with a defined phase 
difference giving rise to so-called metachronal waves. As 
reviewed in the introduction, analytical and numerical in- 
vestigations of beating cilia suggest that hydrodynamic 
interactions between single cilia synchronize their beat- 
ing by inducing such a phase lag and therefore waves to 
occur. For example, in Ref. |17| model cilia are driven 
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pumping 
performance 
per cilium 



a) 




in-line 



Fig. 8. Effective pumping performance ^e// for a two-cilia 
system. One cilium is placed ai x — y — 0, the other cilium 
at different points in the xy plane with coordinates in units 
of L. The effective pumping performance per cilium, ^e//, is 
given in units of the performance ^ for a single, unperturbed 
filament actuated with the same parameters Sp — 3, Bs — 3 
and e — 0.9. The base contour lines run from 0.4 (innermost) 
to 0.9 (outermost) in steps of 0.1. The arrow indicates the 
pumping direction along the y axis. 



by imposing a torque at the base of each ciUum. This 
torque is reversed once the angular amplitude of the cil- 
ium close to the base exceeds a certain threshold value. 
Due to this mechanism, the phase of the driving torque 
and therefore of each cilium can change during the beat- 
ing cycle and a phase lag between neighboring cilia can 
evolve. This cannot happen for the magnetically actuated 
cilia studied in this paper since the oscillating magnetic 
field determines the phase of each cilium. Nevertheless, 
one might ask if in our system the zero-phase-lag state is 
stable. We have never observed an instability towards a 
state with non-zero phase lag. This seems to be reason- 
able: the actuating magnetic field ultimately determines 
the forces which drive the filament against the hydrody- 
namic friction through the liquid. Since these forces are 
much larger than the hydrodynamic interactions between 
the cilia, a noticable phase lag between the cilia cannot oc- 
cur. On the other hand, the artificial cilia system can be 
used to investigate in a systematic way how the pumping 
performance depends on the phase lag between the cilia 
by actuating each cilium with a separate magnetic field. 
In the following, we present initial results of such a study. 

In Fig. [5] we first investigate the pumping performance 
of a two-cilia system with zero phase lag between the cilia. 
One cilium is placed at the center x = y = and the 
other at different points in the xy plane. We introduce an 
effective pumping performance per cilium, ^e//, in units 
of the performance ^ of a single cilium actuated with the 
same magnetic field cycle. If the cilia are sufficiently well 
separated from each other, ^e// approaches the pump- 
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phase shift A(p (in units of 271) 




phase shift A(p (in units of 271 ) 

Fig. 9. a) Average effective pumping performance ^e// for 
multi-cilia systems in the in-line configuration as a function of 
the phase shift between neighboring cilia, ^ef f is given in units 
of the performance ^ for a single, unperturbed filament actu- 
ated with the same parameters Sp = 3, Bs = 3 and e = 0.9. 
b) Average dissipated power Vef / for multi-cilia systems in the 
in-line configuration as a function of the phase shift between 
neighboring cilia, i/^ff is given in units of the dissipated power 
v for a single, unperturbed filament. 



ing performance of a single filament, as expected. If the 
cilia approach each other, £^eff decreases. Due to hydro- 
dynamic interactions between the cilia beating with the 
same phase, their effective friction with the surrounding 
fluid is reduced and therefore they pump less fluid. Figure 
|8] clearly shows that this effect is more pronounced in the 
parallel configuration, where the cilia are placed next to 
each other on the x axis while the strokes are performed 
in the yz plane. On the other hand, in the in-line configu- 
ration, where the cilia sit behind each other on the y axis, 
^eff has almost reached the single-cilium performance ^ 
at a distance L. 

Next we study the pumping performance and the dis- 
sipated energy of several cilia placed on the y axis in 
the in-line configuration. Neighboring cilia are actuated 
by different magnetic fields with a prescribed phase shift 
A(p. As illustrated by the first row in Fig. fTOb). a cil- 
ium lags behind its neighboring cilium to the left by the 
small phase shift A(p. Note that the slow transport stroke 
with a straight cilium goes to the right as can be seen 
from the snapshots at different times in Fig.fTUk). Further- 
more, when time proceeds, the strongly bent cilium in the 
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recovery stroke also moves to the right. So metachronal 
wave propagation and transport stroke occur in the same 
direction which is also termed symplectic metachronism 
in literature [2]. However, note that in biological systems 
the transport stroke is faster than the recovery stroke. 
On the other hand, for phase shifts close to 2tt [see Fig. 
[TOb)] a cilium moves ahead of its neighbor to the left. 
As a result, the metachronal wave propagates opposite to 
the transport stroke and is therefore termed antiplectic. 
Figure [9^) demonstrates that the average effective pump- 
ing performance per cilium, ^e//, falls below the refer- 
ence value of a single cilium for small A(p while it as- 
sumes a maximum value at around Aip = 0.8 x 2tt, i.e., 
in the antiplectic mode. The maximal pumping perfor- 
mance ^e/ / increases with the number of cilia and slightly 
moves to larger Aip. Even at the relatively large distance 
of 1.5L between the cilia, ^e// is increased by more than 
40% relative to a single cilium. Preliminary calculations 
show that this value strongly increases, when the cilia arc 
moved closer together. Hence, our study clearly demon- 
strates that metachronal coordination of ciliary beating 
at the right phase shift significantly enhances the ability 
to transport fluid. An explanation for this behavior can 
be inferred from the snapshots in Fig. 1101 In the optimum 
stroke [see Fig. [TOb)]. the fifth cilium from left performs 
the recovery stroke against the neighboring fourth cilium 
which hinders the fluid flow initiated by recovery stroke 
and therefore increases fluid transport to the right. On the 
other hand, in the metachronal wave with lowest pump- 
ing performance [see Fig.fTUb)]. the cilium performing the 
recovery stroke is further away from the neighboring cilia 
and therefore its fluid transport is hardly hindered. While 
the pumping performance is strongly influenced by the 
phase shift and the number of beating cilia, the average 
dissipated power per cilium, Veff, only exhibits very weak 
dependence on these parameters. As illustrated in Fig. 
ISt) , i^ef f deviates from the single-cilium value by at most 
1 %. This is in contrast to the work of Gueron and Levit- 
Gurevich, who, e.g., found that for a 10-cilia system with a 
distance of 0.7 L between the cilia i^eff decreases by ca. 40 
% relative to the value of a single cilium [5T] . We checked 
that the smaller distance is not the major reason for this 
difference. It might be due to the fact that these authors 
model biological cilia with their internal actuation mech- 
anism, whereas our cilia are actuated by an external field. 



6 Conclusion 

We have performed a thorough theoretical investigation 
of a magnetically powered low Reynolds number pumping 
device. It consists of a superparamagnetic elastic filament 
that is attached to a surface and that can be conveniently 
actuated by an external magnetic field. We have intro- 
duced a reduced pumping performance ^ that quantifies 
the fluid transport and compares it to the value of an ide- 
alized transport stroke of a rigid rod. With the help of ^, 
we have identified an optimum stroke pattern of the mag- 
netically actuated artificial cilium. It consists of a slow 
transport stroke and a fast recovery stroke in contrast 



a) A(p = 0.08*271 
= 0.56 r 
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= 0.66 r 
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A(p = 0.8*271 
= 0.56r 



= o.6ir 
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Fig. 10. Snapshots of the in-line configuration with eight cilia 
for different phase shifts; a) Aip = 0.08 x 2n (minimal pump- 
ing performance) and b) Aip = 0.8 x 2n (optimum pumping 
performance). The times given in units of the time period T 
belong to the field cycle of the rightmost cilium. 



to nature. Our detailed parameter study for the reduced 
magnetic field Bs and the sperm number Sp reveals that 
for sufficiently large Bg the behavior of the artificial cilium 
is mainly governed by the Mason number that compares 
frictional to magnetic forces. Initial studies on multi-cilia 
systems reveal that the pumping performance is very sen- 
sitive to the imposed phase lag between neighboring cilia, 
i.e., to the details of the initiated metachronal waves. Due 
to our results, we expect that metachronal waves in real 
cilia systems also increase the pumping performance for 
fluid transport. 

Experiments using the superparamagnetic elastic fila- 
ment for fluid transport are currently under way |50) . Our 
study will help to identify relevant parameter ranges in 
which the artiflcal cilium can be operated optimally. Initi- 
ating metachronal waves in the multi-cilia system will be 
important for increasing the pumping performance. This 
requires oscillating magnetic fields the phases of which 
should vary from one cilium to the other. The realization 
of such a system is not unrealistic ^51J but certainly poses 
a challenge to experimentalists. It would help to establish 
the artificial cilium as a useful tool for fiuid transport in 
microfluidic devices. 
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A Simulation Parameters 

In table [U we summarize the parameters used for the nu- 
merical simulations. The oscillation frequency to = 2^ /T 
and the magnetic-field strength B were varied to study 
the respective ranges of sperm number Sp and reduced 
magnetic field Bg. The spring constant k, the time step 
At for the Euler integration, and the number Ua of sim- 
ulation cycles were adjusted as necessary. Typical values 
are shown. 



B Mobilities 

Explicit expressions for the stresslet and sourcelet contri- 
cutions of Blake's tensor can be calculated from Ref. [53] 
or taken directly from Ref. [52| . Starting with these, we 
have calculated the self- and cross-mobilities as outlined 
in the main text. 



B.l Self mobilities 

Due to the axial symmetry of the sphere with its image 
system all but the diagonal elements vanish and the xx 
and the yy elements are identical: 



Table 1. Parameters used in the numerical studies. 



SfJ-self 




(32) 



where u is given by (with a being the radius of the bead 
and z the distance to the plane) 



3 
16 



1 /a 
3 Vz 



(33) 



Once the expression above is included into formula ([T^ 
for /iji, the self-mobility agrees with the expression given 
in Ref. HTl. 



B.2 Cross mobilities 

Symmetry dictates that the xy and yx elements are iden- 
tical (the boundary wall is in the xy plane). Furthermore, 
the XX and yy elements are formally equivalent when 
is interchanged with Ry. The same holds for the xz and 
yz or zx and zy elements, respectively. 

The positions of the beads are ri = (xi,yi, zi) and 
T"2 = {x2, 2/2 7 ^2) and the image system is located at ¥2 = 
(2^2, 2/2, — •22)- Furthermore, for convenience we define s — 
\ri - r2\ and R^ = {xi - X2), Ry = {yi - 2/2) and R^ = 
(zi -t- Z2). Let a and f3 refer to the x and y but not the 
z coordinate. In this notation, we obtain the following 
components for the cross mobilities: 



Parameter 


Simulation value 


N 


20 


a [/xm] 


0.5 


lo 


3a 


— > L [/im] 


28.5 


X 


0.993 


7] [Ns/m^] 


10"^ 


k [N/m] 


1.5 • 10"^ 


A [Nm] 


4.5 ■ 10"^^ 


to [2tv/s] 


0.071 . . . 140.72 


^ Sp 


1.5 ...10 


B [T] 


. . .0.085 




...7 


V^max [ °] 


60 


At [s] 


0(10-'^) 


ris 


f» 5 
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